{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "2a0b97d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.interpolate import interp2d, interpn\n",
    "from numpy.matlib import repmat\n",
    "import math"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00f6be88",
   "metadata": {},
   "source": [
    "# Entry/Exit with Homogeneous firms\n",
    "## and translog demand\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "244c912c",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sequence_jacobian import solved, simple, create_model\n",
    "\n",
    "@simple\n",
    "def household(w, psi, nu):\n",
    "    L = (w/psi)**nu\n",
    "    return L\n",
    "\n",
    "@solved(unknowns={'v': 1.0, 'N': 1.0}, targets=['valfunc', 'accumulation'])\n",
    "def firms(C, sigma, v, beta, delta, Ne, N, L):\n",
    "    mu      = 1 + 1/(sigma*N)\n",
    "#    y       = p**(-sigma)*C\n",
    "    ell     = L/N            # imposing labor market clearing\n",
    "    pi      = (1-1/mu)*(C/N)\n",
    "    valfunc = v - (pi + beta*(1-delta)*v(1))\n",
    "    accumulation =  N - (1-delta)*N(-1) - Ne\n",
    "    return mu, ell, pi, valfunc, accumulation\n",
    "    \n",
    "@simple \n",
    "def mkt_clearing(mu, N, v, fE, w, Z, C, L, pi, Ne, Nbar):\n",
    "    p       = mu*w/Z\n",
    "    rho     = math.e**(-0.5*(Nbar - N)/(2*Nbar*N))\n",
    "    goods_mkt = p - rho                   # as in BGM\n",
    "    free_entry = v - fE\n",
    "    budget     = C - w*L - N*pi\n",
    "    return goods_mkt, free_entry, budget, rho\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0a70a25",
   "metadata": {},
   "source": [
    "### Steady State"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "d3fbc144",
   "metadata": {},
   "outputs": [],
   "source": [
    "ces = create_model([household, mkt_clearing, firms], name=\"CES Model\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "bb885ca8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.1839464882943145"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "calibration = {'beta':    0.96,\n",
    "               'psi':     1,\n",
    "               'fE':      .6339,\n",
    "                'nu': 2.0,\n",
    "                'delta':  0.11,\n",
    "                'sigma':  2.6,\n",
    "                'w': .75,\n",
    "                'Z': 1.01,\n",
    "                'Ne': .23, \n",
    "                'C': 1, \n",
    "                'Nbar': 2}\n",
    "\n",
    "ss = ces.steady_state(calibration)\n",
    "ss['mu']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "283c6fca",
   "metadata": {},
   "outputs": [],
   "source": [
    "w0    = .75\n",
    "Nbar0 = 1/(calibration['sigma']*(1/w0-1))\n",
    "Ne0   = calibration['delta']*Nbar0\n",
    "\n",
    "calibration['Ne'] = Ne0\n",
    "\n",
    "ss=ces.solve_steady_state(calibration, \n",
    "      unknowns={'Nbar': Nbar0, 'w': w0, 'Ne': Ne0}, \n",
    "      targets={'rho': 1.0, 'free_entry': 0.0, 'goods_mkt': 0.0})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "2e4bcb7f",
   "metadata": {},
   "outputs": [],
   "source": [
    "ss=ces.solve_steady_state(calibration, \n",
    "      unknowns={'Nbar': ss['Nbar'], 'w': ss['w'], 'Ne': ss['Ne'], 'C' : 1.0}, \n",
    "      targets={'rho': 1.0, 'free_entry': 0.0, 'goods_mkt': 0.0, 'budget' :0.0})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "051e3acb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.2636807129968646"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ss['mu']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "8240f673",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.458640566630893"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ss['N']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2e52b213",
   "metadata": {},
   "source": [
    "# Shock to sunk cost of entry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "1d08066a",
   "metadata": {},
   "outputs": [],
   "source": [
    "T = 1000\n",
    "dZ = .68*.02 * .685 ** np.arange(T)\n",
    "irf = {}\n",
    "irf = ces.solve_impulse_linear(ss, ['w', 'Ne', 'C'], ['free_entry', 'goods_mkt', 'budget'], {'fE': dZ})\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "c00e2675",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAEYCAYAAABRMYxdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABO/klEQVR4nO3deXxddZ34/9f73pu9Sbpka5O06U73hUBZKpRNFhVcEUQGGJHRUUfHBdcZnfk6o6P+3JlBEFARxYoiiAiURQrI0hYaugJd0iZNmq3Nvifv3x/npFzSLDfLzbn35P18PO4j9yw5533Szzl9n8/5nM9HVBVjjDHGGD8JeB2AMcYYY8x4swTHGGOMMb5jCY4xxhhjfMcSHGOMMcb4jiU4xhhjjPEdS3CMMcYY4zuW4ESRiHxcRKpEpFlEZngUQ5GIqIiEBln+DRH5tft9thtrcJhtXi8iz0YjXmMiMVy5NsaYuExwRKRURDpFJKvf/O3uRa/Io9DCY0kAvg+8XVWnqGrdOG03asmFqh52Y+2Jxvajwf33XuB1HGZw8XC+GjPR3JvJvk+viLSFTV/j3nx29VvvZvd3/yYi7e68WhH5o4jM9PqYYk1cJjiug8DVfRMisgJI8S6ck+QCycAurwMxJgaM2/lqtTbGD9ybySmqOgU4DLwrbN497mq/C19PVb8TtolPur+7CJgK/GBijyD2xXOCczfwD2HT1wG/Cl9BRN4hIq+ISKOIlInIN8KWJYvIr0WkTkTqRWSLiOS6y64XkQMi0iQiB0XkmoECEJEkEfmhiFS4nx+68xYBr7mr1YvIk4P8/hki8nd3/yUisiFs2UkxiMgS4FbgTDdzrx/uOMP8oxtjpYh8bpB43lLtP9zfQUS+JyLH3WWXhs3/m4h80z22ZhH5s4jMEJF73Bi3hN+1i8gpIrJJRI6JyGsicmXYsl+IyC0i8hc3jhdFZL67bLO7Wom7nw+KSJaIPOT+TY+JyDMiEs/l3C+GPF+HOVf7yuVHROQwcNL5JCLvc2uKlrtl5pthyzaISHnYdKmIfFlEdrvl9y4RSR7n4zVmQqjqMeAPwHKvY4k18XzhfwHIEJEl4rQZ+SDw637rtOBcVKcC7wA+LiLvdpddB2QChcAM4GNAm4ikAT8GLlXVdOAsYPsgMXwVOANYDawCTge+pqqvA8vcdaaq6vn9f1FE8oG/AN8EpgOfB/4gItmDxaCqe9w4n3ez+akRHGef84CFwNuBL4nIhYMcU198w/0d1uEkcVnAd4A7RETCll8FXAvkA/OB54G73GPdA3w9bD+bgN8AOTh3+f8rIsvCtnU18B/ANGAf8F8AqnqOu3yV+/f4HfA5oBzIxqlF+wpg45F4b7jzNZIyfC6wBLg4fKaI3AD8D3Chqu6MMJ5r3O3Mx7kD/tpIDsaYWCHOo9/3Aa94HUusiecEB968K7wI2AscCV+oqn9T1R2q2quqrwK/xblIAnThJDYLVLVHVbepaqO7rBdYLiIpqlqpqoM9ZroG+E9VrVbVGpz/hK+NMPYPAw+r6sNufJuArcBlI4xhuOPs8x+q2qKqO3ASjatP2tDJhorhkKre7rbX+SUwEyeh6HOXqu5X1Qbgr8B+VX1cVbuB3wNr3PXeCZSq6l2q2q2qL+Pcjbw/bFt/VNWX3N+9ByehHEyXG8scVe1S1WfUBlyLFYOerxGW4W+4ZbgtbN5ngC8AG1R13whi+amqlrl3v/9FZOeDMRPtSrc2uu8zK2zZj91a/BKgEvisJxHGMD8kOB8Crqff4ykAEVknIk+JSI2INODUfmSF/e6jwL3uo5vviEiCqrbg3F1+DKh0H42cMsj+ZwGHwqYPufMiMQf4QHjhBdYDM0cYw3DH2adsJHFGEMPRsHVb3a9TwpZXhX1vG2C6b905wLp+f4drgLyB9gW09ttPf9/FqeV5zH289qUh1jUTa9DzdRRluM8XgFtUtXyAZUMZ0flgjEc2qurUsE9F2LJ/ceflq+o17k22CRPXCY6qHsJpvHgZ8McBVvkN8CBQqKqZOO1XxP3dLlX9D1VdivP45Z24bQRU9VFVvQinJmAvcPsgIVTg/AfdZ7Y7LxJlwN39Cm+aqn57mBgGqo0Y9DjDFI40zhH8HcaiDHi6399hiqp+fDQbU9UmVf2cqs4D3gV8VkQuGNeIzagMc75GUoYHKvtvB74mIu8Lm9cCpIZN53GyEZ8Pxpj4EtcJjusjwPlujUN/6cAxVW0XkdNx7h4BEJHzRGSF2x6gEefRRo+I5IrI5W7bkA6gGRjstenf4lxcs93noP/Oye2ABvNr4F0icrGIBMVp9LxBRAqGiaEKKBCRxEiOM8y/iUiq27blBuB3QwU3wr/DWDwELBKRa0Ukwf2cJk6D6khUAfP6JkTknSKywG0P1OjGHDevvU8Cg52vkZThgewCLgFuEZHL3XnbgctEZLqI5OE8xurvE+65Nh2nndaQ54MxJv7EfYLjtvPYOsjifwb+U0SacJKPjWHL8oD7cP4T3AM8jZN0BHAaqlYAx3DaAfzzINv/Jk67mVeBHcDL7rxI4i4DrsC5uNbg1GR8wd3/UDE8iXNRPyoitREcZ5+ncR7dPAF8T1UfGybEkfwdRk1Vm3Duwq9y93UUp8FoUoSb+AbwS/fx1pU4Dakfx0nIngf+V1X/Ns5hm1Ea4nyNpAwPts0SnBrY28V5m+9unHYJpcBjDJy8/MZddsD9RHTeGmPih1j7S2PMZCIipcCNqvq417EYY6In7mtwjDHGGGP6GzbBEZGz3XYYiMiHReT7IjJnuN8zxhhjjPHKsI+oRORVnE7sVuI8274DeK+q9u+jwhhjjDEmJkTyiKrb7SjtCuBHqvojnDcejDHGGGNiUiSD1jWJyJdxet49x32tOiG6YQ0sKytLi4qKvNi18Ylt27bVqmq213EMx8q6GSsr62ayGKysR5LgfBCnT4qPqOpREZmN01vshCsqKmLr1sHeCDdmeCJyaPi1RrS9O3FeUa5W1ZMGuxNnANUHcDq4A2fYif8cbrtW1s1YjXdZjxYr62asBivrkTyiasJ5NPWMOKNkr8bp4G48grpEnNGj91mX+iZO/QKno7mhPKOqq93PsMmNMbFouOu1OH7sLn9VRNZ6EacxfSKpwdkMvE1EpuF0ErcVp1bnmrHs2H3UdQvOwHvlwBYReVBVd49lu36iqqg6/dOrqvsT1O2xfrgujCZrF0cikJwQnJB9qepmESmakJ0BzR3dJIcChILWw4OZOBFery/F6WhzIbAO+D/3pwnTd13vdS/Qdl2PzGiu65EkOKKqrSLyEeAnqvodEdk+mgD7OR3Yp6oHAETkXpyGzIMmOHV1dWzfvp3Vq1fT09PD3Xffzdq1a1m5ciVdXV3cc889FBcXs3z5ctrb27n33ntZt24dS5YsobW1lY0bN3LmmWeyePFimpubue+++1i/fj0LFiygoaGB+++/n3POOYd58+Zx/PhxHnjgATZs2EBRURG1tbU89NBDXHDBBRQWFlJdXc3DDz/MhRdeSFJmNrv3H2bbc0+RsaCYlmA6x2uO0lNWQt30ZTSQgrTUkdu0j71JC2nUZNI765nfc4htzKepN4kZWs9SOcJz3XNp1iRmSj3Lg5U83TmfNhIoDNSzLHSUpzrn00ECcwLHWRKq4onOBXQRYm7wGIuD1WzqXEgPQeYF61gUrOHRzkUoARYEa1kQrOWRTme8zEXBGuYGj/Fo52IATglWUxisZ1PnIgCWBquYGWzkic6FACwPHSU70MxTnQsAWBGqZLq08nTXfABWhSrIlHY2dzmjJqwJHSFNOnm2ay4Ap4bKSZJu/t5VBMBpoTKC0ssLXU6PA6cnHAbgpa7ZAJyRcIgeDbCl2xky6KyEUjo0xLbuAgDWJxykRRN5pTsfgHMSDtCgyZR0z2JmZjKfW3icvLw81q9fD8DGjRspKCgYXUkduzNFpASnp+bPDzYyvIjcBNwEMHv27JOWP7Gnio/8cisPfvJsVhZMjWK4xpwkkuv1FcCv3JdSXhCRqSIyU1Urh9pwLF7XN5x/IT0pU3njYBk7XtxM8ty1tATTaayrQip2UJO5hHpSCbXWkdeyn72JC2noTSa9+zgLew6zVefR2JtEFvUskwqe6ZpLsyYySxpYEbLr+miu66fOmcbV2RUjuq5HlOCIyJk4NTYfceeNx+1xPm8d0becAbL98It+fn7+OOx2bKoa29n6Sjnb9pbSUdHAXT9/kUPtyUyXVk5PaOKl8v0c01RmJ7WzJthNZX0bkprI9FCQxFCA+VlTCKZmkNSuJNUkctHsHAIp6QSaE9HKOq6eP4dQchq99ZV0VjTw0cVzCSam0n2snPYjjXxi6XyCiSl01B6m7UgTn16+kEAokY7qQ7RVNPO5FYuRYIj2qgO0VbZw86pTkECAtsr9tB1t5UtrnBOhtSJIR3UHX1rtTh8ROuq6+dJKd7ocOo73ctoKZ7rlcC9djcK65e70oW66mo9z5jJnurm0k57WRs5a6k4f7KCno4X1pzjTTQfa0K4OzlnsTu9rQXt72LDInX6jGYDzFzrTja83IoEgFyxwp1+rRxKSuGieM92w9xjBpDQunutO764lmJrBpUWnkJYUgvItUfn3H4WXgTmq2iwilwF/wrnDPYmq3gbcBlBcXHzSfVrWFGf0iqrGjmjFasxgIrleD7ROPnBSghNL1/VeVXYcqefBfR3s3FdG0tF6brn179T0prnX9WZeqiilKZDG7OQOVtLN8dYuAqmQnpJAcmeARbnpBFMySGqHhOqjXDwnj2BKOtKUTG/lMT68wLmu99RX0nmkgZsWzyOYlEJXnV3XI72u52Uk071/ZGPiRtIPzjnA54HnVPV/RGQe8BlV/ZcR7enk7X4AuFhVb3SnrwVOV9VPDfY7xcXF6kVjtAM1zfy5pJK/7Kjg9SrnH2xKUojl+RksyJnC/OwpFE5LJS8zmZyMJKalJpJgjxBikohsU9Xicd5mEfDQQI2MB1i3FChW1dqh1huorB9taOeMbz3Bf71nOdess742zdDGs6xHcr0Wkb8A31LVZ93pJ4CbVXXbUNv24rre1dPLY7uqeHhnJX/bW01LpzMe77ysNJbMymB+VhpFWWnMmppCXkYy2elJpCYGccbwNbFmsLI+bA2Oqm7GaYfTN30AGFNy4yoHCsOmC3Cq8GPGziMN/OTJN3h0VxUicFrRdL7+rqWsmzuDxXnpBANW2M3Q3NGsq1RV3VGyA0DdaLaVNSUREai2Ghwz8SK5Xsf8Nb2lo5t7t5RxxzMHqGhoJ2tKIpevzufCJTmsnT2NaWmJXodoxlEkj6iiZQuwUETmAkdwRpP+kIfxnNDQ2sXXHtjJn0sqSE8O8S8XLOSadbPJzUj2OjQTY0Tkt8AGIEtEyoGv4/YTpaq3Au8HPi4i3UAbcJWOcoTbUDDAjLREqpvaxyV2Y0Ygkuv1g8An3fY564CG4drfTKSn9lbzxT+8SnVTB6fPnc5/vWcF5yzKthtVH/MswVHVbhH5JPAoTpueOwdrfDmRXjhQx7/+bjs1TR18+oKFfORtc8lI9qRfQxMHVPXqYZb/FPjpeO0vJz3ZanDMhBvsei0iH3OX3wo8DFwG7ANagRu8ijdcc0c333xoN/duKWNxbjr/9+G1nDpnutdhmQngZQ0OqvowzkkRE+5+vpR/f3AXRTPSuP+fz2ZFQabXIRnzFjkZSVQ3WYJjJt5A12s3sen7rsAnJjquoTS0dXHtHS+y80gDHzt3Pv960UKSQhPThYTx3rAJjlsl+SmgKHx9Vb08emFNvI1byvi3B3Zx4ZIcfnTVGudNHGNiTE56ErsrGr0Ow5iY15fc7Kls5GfXFnPR0lyvQzITLJL/xf+EM4L4n4HeqEbjkQdLKvjiH1/lbQuzuOWatZbhm5iVm5FMbXMHPb1qbQeMGURzR/eJ5Ob/rjmVCy25mZQiSXDaVfXHUY/EIzuPNPDZ323ntKLp3HZtsSU3JqblpCfRq1DX0kFOujV6N2Yg33xoNzuONHD7tcWW3ExikSQ4PxKRrwOPASce/qvqy1GLaoJ0dvfy+d+XMD0tkduvLSYl0ZIbE9uy3aSmutESHGMG8uTeKu7dUsY/nTvPkptJLpIEZwVwLXA+bz6iUnc6rv3f3/az92gTt/9DMZmp9qaUiX25GU5vxs6r4tYI3phwx1s6+eIfdrA4N53PXrTI63CMxyJJcN4DzFPVzmgHM5H2Hm3kp0+9weWrZlnjMxM3cjLerMExxrzV1x/cRX1rJ7+44TRrbmCIZDyBEmBqlOOYcP/2p51kJCfwjcuXeR2KMRHLtvGojBnQ3qONPFhSwcfOnc+yWVa7aSKrwckF9orIFt7aBiduXxPfdugYW0qP8/V3LWW6dc1t4khiKMB0683YmJP87OkDpCYG+cj6uV6HYmJEJAnO16MexQS7ffNBMpJDXFlcOPzKxsSYnHTr7M+YcOXHW3mwpILrzypiaqrdtBrHsI+oVPVpYC+Q7n72uPPi0qG6Fh7dfZRrzphjnfmZuJSdnkR1o9XgGNPn588cRMBqb8xbDJvgiMiVwEvAB4ArgRdF5P3RDixa7nqulFBAuP6sIq9DMWZUcjOSrQbHGNexlk7u3XKYK1bnM2tqitfhmBgSSRXGV4HTVLUaQESygceB+6IZWDQ0tHaxcWsZ71o1y0YGN3ErJz2JmqYOenuVgPVmbCa5X79wiPauXj527jyvQzExJpK3qAJ9yY2rLsLfizm/31ZGa2cPN663E8HEr5z0JLp7leOtvuq5wZhR+XNJBevmTmdhbrrXoZgYE0mi8oiIPCoi14vI9cBfiKERwEfikZ1HWTozg6WzMrwOxZhR66t9tFfFzWR3oKaZN6qbuWR5ntehmBg0ZIIjIgL8GPgZsBJYBdymql+cgNjGVV1zB9sOH7dO/Uzcy3lLb8bGTF6P7a4C4O3LLMExJxuyDY6qqoj8SVVPBf44XjsVke8C7wI6gf3ADapaP17bH8gTe6tRxRIcE/f6xqCyhsZmsnt011FW5GeSb42LzQAieUT1goicNs773QQsV9WVwOvAl8d5+yfvcHcVszKTWWaPp0ycy053a3DsVXEziVU1tvPK4XouXmY3rWZgkSQ45wHPi8h+EXlVRHaIyKtj2amqPqaq3e7kC0DBWLY3nLbOHp55o4YLl+biPHUzJn4lJwTJTEmwGhwzqfU9nrrYHk+ZQQz6iEpE5qrqQeDSKMfwj8DvhojjJuAmgNmzZ49qB8/uq6W9q9ceTxnfyElPsgE3zaT22K6jzMtKY0HOFK9DMTFqqBqcvn5u7lTVQ/0/w21YRB4XkZ0DfK4IW+erQDdwz2DbUdXbVLVYVYuzs7MjPa632LT7KOlJIdbNnTGq3zdmMCJyp4hUi8jOQZaLiPxYRPa5NaBrx2O/uRnJVFkjYzNJNbR28fz+Ot6+LM9q5c2ghmpkHBCRrwOLROSz/Req6veH2rCqXjjUchG5DngncIGqaiTBjkZPr/LEnmo2nJJDYiguu+8xse0XwE+BXw2y/FJgoftZB/yf+3NMctKTePFgy1g3Y0xc+vv+Wrp7lYuW5ngdiolhQ/2PfxXQjpMEpQ/wGTURuQT4InC5qraOZVvD2V3RSF1LJxecYieCGX+quhk4NsQqVwC/UscLwFQRmTnW/RZMS6GyoY3O7t6xbsqYuLO9rJ7EYIDl+Zleh2Ji2KA1OKr6GvA/IvKqqv51nPf7UyAJ2ORWL76gqh8b530AsL28HoBT50yLxuaNGU4+UBY2Xe7Oq+y/4kjam82ZkUavQtnxVuZnWxsEM7lsL6tn6awMkkJBr0MxMWzYsaiikNygqgvGe5uDKSmrZ0ZaIgXTrJ8E44mBGggM+EhWVW8DbgMoLi4e8rFtUVYaAKW1LZbgmEmlp1fZcaSBK4sLvQ7FxLhIBtuMayVl9awqnGoN0YxXyoHwK3EBUDHWjc51E5yDtdYOx0SfiEzHedu1CCgFrlTV4wOsVwo0AT1At6oWj3csr1c10drZw+rCqeO9aeMzvm5129Texb6aZlYVTPU6FDN5PQj8g/s21RlAg6qe9HhqpKalJpCRHKK0zhIcMyG+BDyhqguBJ9zpwZynqqujkdyAc9MKWIJjhhVRDY6InIWTuZ9YX1UHe2skZuw40oAqrCq0hmgmOkTkt8AGIEtEyoGvAwkAqnorzsC0lwH7gFbghnHaL3Oz0iitjWobfWP6XIFTzgF+CfwN50WRCbe9rJ6pqQnMmZHqxe5NHBk2wRGRu4H5wHacakdw2hDEfIJTUtYAYDU4JmpU9ephlivwiWjsuygrja2lJz0lMCYacvtqHlW1UkQGey1VgcdERIGfue3KTjKWDly3l9WzqsCaHZjhRVKDUwwsjWZfNdFSUlbPnBmpTEtL9DoUY8Zd0Yw0HiypoL2rh+QEe5vEjI2IPA4MNO7BV0ewmbNVtcJNgDaJyF63K4W3GEmD+nAtHd28XtVkwzOYiESS4OzEKfRjbjcw0UrK6zmtaLrXYRgTFXOz0lCFsmOtLMwdU9dUxgzZOauIVInITLf2ZiZQPcg2Ktyf1SJyP3A6cFKCM1o7jjTQq9b+xkQmkkbGWcBuEXlURB7s+0Q7sLGqamynsqGdVXYiGJ8qsjepzMR5ELjO/X4d8ED/FUQkTUTS+74Db8e5QR43290GxnZdN5GIpAbnG9EOIhrebGlvDYyNP82d4faFY29Smej7NrBRRD4CHAY+ACAis4Cfq+plQC5wv9s2JgT8RlUfGc8gSsrqmT09lenW7MBEIJKO/p4WkVzgNHfWS6o6YPVkLCkprycYEJbNsgTH+FNmagLTUhM4aG9SmShT1TrgggHmV+C8JYiqHgBWRTOOkrJ6iq3ZgYnQsI+oRORK4CWcjP1K4EUReX+0AxurkrIGTslLt8aXxteKstIotUdUZhJoau+ioqGdJTMzvA7FxIlIHlF9FTitr9ZGRLKBx4H7ohnYWO092sj5NsCm8bm5M9J4/kCd12EYE3WH6pyayrlZ1v+NiUwkjYwD/R5J1UX4e55pau+itrmTuVk2Ro/xt6KsNCob2mnr7Bl+ZWPiWF9j+r7G9cYMJ5IanEdE5FHgt+70B3F6Z41Zfb27WqZv/K7vYn/oWAun5FnVvfGvvkexc6ZbgmMiE0kj4y+IyPuAs3FGRr5NVe+PemRjcLDOMn0zOZx4k6rWEhzjb6V1reRlJJOSaO0qTWQiGotKVf8A/CHKsYwby/TNZFHk1lLam1TG70rrWk6Ud2MiMWhbGhF51v3ZJCKNYZ8mEWkcj52LyOdFREUkazy216e0toWZmZbpG/9LT04gOz2JfdXNXodiTFSV1rZQNMNuWk3kBq3BUdX17s+o9AEvIoXARTidRo2rg3V2IpjJY0V+JjuO1HsdhjFR09jeRV1LpzU7MCMSST84d0cybxR+ANyMM/rsuCqtbbETwUwaqwqm8kZ1M80d3V6HYkxUHHIfwdqNqxmJSF73XhY+ISIh4NSx7FRELgeOqGpJBOveJCJbRWRrTU3NsNtuaO3ieGuXvUFlJo1VhZmowo7yBq9DMSYq3nxxxK7rJnJDtcH5sog0ASvD298AVQww0NoAv/+4iOwc4HMFTueB/x5JgKp6m6oWq2pxdnb2sOufOBEs0zeTxKqCqYAzPIkxfnTIXhwxozBUG5xvAd8SkW+p6pdHumFVvXCg+SKyApgLlLiDshUAL4vI6ap6dKT76a/vDaq59ojKTBLT0hKZMyP1xACzxvjNwTp7ccSMXCT94HxZRKYBC4HksPmbR7NDVd0BnBhDQURKgWJVrR3N9vo7WNuCCBROt6pMM3msKpjK1tJjXodhTFSU1rYwZ4Zd083IRNLI+EZgM/Ao8B/uz29EN6zRK61rYVZmig2yaSaVVYVTqWhop7qx3etQjBl3h+parVbejFgkjYw/DZwGHFLV84A1wPCtfSOkqkXjVXsDTqZvJ4KZbFYXZgJQYg2Njc+ceEXc2lWaEYokwWlX1XYAEUlS1b3A4uiGNXqlda3W0t5MKBG5REReE5F9IvKlAZZvEJEGEdnufiJqYD8Sy2ZlEgyItcMxvnOiZ3pLcMwIRTJUQ7mITAX+BGwSkeNARTSDGq3jLZ00tHVZpm8mjIgEgVtwOq0sB7aIyIOqurvfqs+o6jujFUdyQpBT8tLtTSrjO6V1fYMn23XdjEwkjYzf4379hog8BWQCj0Q1qlHqe0XcTgQzgU4H9qnqAQARuRe4Auif4ETdqsKpPFRSQW+vEgjIRO/emKh4swbHaubNyAzVD06G+3N63wfYATwLTJmg+Eak70SwXozNBMoHysKmy915/Z0pIiUi8lcRWTbA8hF3atnf6oKpNLZ3U+om+sb4Qan7iri9OGJGaqganN8A7wS24QynIP1+zot6dCNUWtdKQKBwmmX6ZsIMVFXSf/iRl4E5qtosIpfhPO5deNIvqd4G3AZQXFw84iFM1s6ZBsBz++uYlx2T9yDGjFh1YwczM5OHX9GYfgatwelrL6Cqc1V1Xv+fExdi5I42tJGdnkRiKJK208aMi3KgMGy6gH5t1FS1UVWb3e8PAwkikjXegczPTqNoRiqbdleN96aN8Ux1UzvZ6Uleh2HiUCT94DwgIleLSMxXi1Q1dpCTbpm+mVBbgIUiMldEEoGrgAfDVxCRPHG77RaR03HOu7rxDkREuGhpLs/vr6WpvWu8N2+MJ2qaOizBMaMSSVXH94G3AXtE5Pci8n4Ricksorqpg9wMOxHMxFHVbuCTOB1g7gE2quouEfmYiHzMXe39wE4RKQF+DFylqiN+BBWJi5bm0dWjbH593LqWMsYznd29HG/tshtXMyqRvEX1NPC0+zrs+cBHgTuBjCjHNmI1Te2sLpzqdRhmknEfOz3cb96tYd9/Cvx0ImI5dc40pqUmsGn3Ud6xcuZE7NKYqKlt7gCwGhwzKhE1VhGRFOB9wMdwejX+ZTSDGo2unl5qmzvJsRPBTGLBgHD+Kbk8ubearp5er8MxPiEiHxCRXSLSKyLFQ6w3ZKeXI1XT5CQ4dl03oxFJG5zf4VS9n4/Todl8Vf1UtAMbqb5MPzfDqjLN5HbR0lwa27vZYoNvmvGzE3gvzriEAwrr9PJSYClwtYgsHctOq5usBseMXiQ9Gd8FfEhVe6IdzFhUN1qmbwzAOYuySAoF2LS7irPmj/vLWmYSUtU94DRkH8K4d3pZYwmOGYNIHlFtBr4sIrcBiMhCEYlal/OjVeWOopxjjYzNJJeaGGL9giwe21VFb29U2jIbM5BIO72MWF+CkzXFrutm5CJJcO4COoGz3Oly4JtRi2iU+qoy7RGVMXDFmnyO1Lfx5N5qr0MxcUJEHheRnQN8roh0EwPMGzDDjrTX7uqmdqanJZIQtL7NzMhF8ohqvqp+UESuBlDVNhmmntIL1U0diMCMtESvQzHGc5cuz2NWZjK3P3OAC5fmeh2OiQOqeuEYNzFsp5dh+4qo1+6apg6yrfbGjFIkaXGn+xaVAojIfKBjrDsWkU+5re13ich3xrq96sZ2ZqQlEbJM3xgSggFuOHsuLx48xo7yBq/DMZPDsJ1ejlR1U4c1OzCjFkk28HWc0cMLReQe4Ang5rHsVETOw2l8tlJVlwHfG8v2wDr5M6a/D55eyJSkELc/c8DrUEycE5H3iEg5cCbwFxF51J0/S0QehsE7vRzLfq0Gx4xFJB39bRKRl4EzcJ6xflpVx9pN6seBb6tqh7uPMTcUqG5qtzeojAmTkZzAVacVctffS/nipaeQPzXF65BMnFLV+4H7B5hfAVwWNn1Sp5dj2Cc1zR1k242rGaVBa3BEZG3fB5gDVOI8T53tzhuLRcDbRORFEXlaRE4bIo6IGqPZOFTGnOyG9XMBuH2z1eKY+NLY1k1nd6/V4JhRG6oG5/9zfyYDxUAJTg3OSuBFYP1QGxaRx4G8ARZ91d3vNJxaodOAjSIyb6DxeSJpjNbTq9Q12yMqY/rLn5rCB08r5FfPl/LuNfk2lImJGzXNTtcf1geOGa1Ba3BU9TxVPQ84BKxV1WJVPRVYA+wbbsOqeqGqLh/g8wBOa/s/quMloBcYdY9kdc0d9Cpk2yvixpzkS5eeQk56MjffV0JHd0z312nMCdUnhmmw67oZnUgaGZ+iqjv6JlR1J7B6jPv9E87QD4jIIiARGHW7nirrxdiYQWUkJ/Df713O61XN3PLksPcmxsQE68XYjFUkCc4eEfm5iGwQkXNF5HacFvJjcScwT0R2AvcC1w30eCpS1U1OVaZ18mfMwM4/JZf3rsnnf/+2n5cPH/c6HGOGZQmOGatIEpwbgF3Ap4HP4IwrcsNYdqqqnar6YfeR1VpVfXIs26u2EWeNGda/v2sps6amcP2dL7GrwvrGMbGtuqmDpFCAjORI+qM15mTDJjiq2q6qP1DV97ifH6hq+0QEF6m+cahsvBJjBjc1NZF7blzHlKQQ197xEvuqm7wOyZhB1TR1kJ2eNNwAn8YMyhfd/lY3dTAjLZHEkC8Ox5ioKZyeyq9vXEdAhKtue5Hn9o21SytjoqOmqcNq5c2Y+CIjqG7ssOe0xkRoXvYUfvvRdWSmhPjwHS/yrYf30Nnd63VYxrxFdVO7XdfNmAyb4IjIByKZ56XqpnZyrIGxMRFbmJvOQ596G1efPpufbT7AxT/czH3byunusUTHxIa+R1TGjFYkNThfjnCeZ6obO8i1E8GYEUlJDPLf71nBXdefRnJCkM//voTz/7+nueWpfRysbfE6PDOJdXb3cry1y/rAMWMyaPN0EbkUZ4yRfBH5cdiiDKA72oFFqrfXGa/ERpw1XhGRS4AfAUHg56r67X7LxV1+GdAKXK+qL094oIM475QcNizO5vE91dy2eT/fffQ1vvvoa5ySl866udM5tWg6qwoyKZiWSjBgDT5N9NU22yviZuyGev+uAtgKXA5sC5vfBPxrNIMaibqWTnp61TJ94wkRCQK3ABfh9NC9RUQeVNXdYatdCix0P+uA/3N/xgwR4aKluVy0NJeK+jYe3lHJk3ur+f22cn75/CEAEkMB5s5Io2BaCnmZyeRmJDMtLZGpKQlkpiSQlhRiSlKIlIQgyQkBkkJBEkMBEoJCKOiL5n5mgtRY1x9mHAya4KhqCVAiIvcDLaraAycu6DFT6vo6+bMTwXjkdGCfqh4AEJF7gStw+ovqcwXwK7czyxdEZKqIzFTVysE2WldXx/bt21m9ejU9PT3cfffdrF27lpUrV9LV1cU999xDcXExy5cvp729nXvvvZd169axZMkSWltb2bhxI2eeeSaLFy+mubmZ++67j/Xr17NgwQIaGhq4//77Oeecc5g3bx7Hjx/ngQceYMOGDRQVFZHY3UJo/2b+55ILmDnrdF7YdZDnn36cztzlHGxLpL62htRDr/HX9gKOaSpZ0kJxQhkvdM2hXlPICTSzNlTO37uKaNRkcgNNrE04wgvdc2kLpDBTGlkWOMJWWUhnIJlcjrOg9wjbg4vpDCSS23uMub1H2B5aQncgkZyeWmb3VLA9YSk9gRB5PTUU9FTySuIyeiVIXnc1+T1HeSVxBSoBZnYfZWZPFS8nrUIEZnVXktNdw/bklQDkd1WQ1VNHSfIKAAq6jjCt9zg7kpYDUNhVRmZvEzuTlgIwp+swU3pb2JW0BICirkOk9raxO+kUAOZ2lpKsHexJWgzAvM6DJGgXryUtAmBB5wEC9PB64kIAFnbuB+CNxPkALOp8g16C7EucB8DijtfpkgQOJDoDpS7peI12SeJgYhEASzv20hpIoTRhDkUz0nhXxmHy8vJYv94ZHnDjxo0UFBSMrBTHmGrr5M+Mg0h6UHoMuBBodqdT3HlnRSuokUhOCPLu1bOYnzPF61DM5JQPlIVNl3Ny7cxA6+QDb0lwROQm4CaA/Pz8cQ90NELBAIvy0tmXnsRFZxeRn5/P0aNHeeSRam6+8GySMqezv7SMrc818q41yyBlKnVVFVTuOc6qRfOQlAxa6ippPXSMD86dTW9iGt31R9HKWi6enUtvQiraCMHqGtYXZtEbTCbQ1E1CbTWnF0xHQ0mEGjtIOBbi1IJpEEggobGNxOMhTi2cBoEQCQ0tJNaHWDt7GkiAxPomEhtCrJ09FYDE+gYSm0KscQcaTTp+nISWBNYU9E0fI9TS/Ob0sVpCbe2syXen62oIdXSyZpYznVxbRaCzO2w6mUCXsmamO12TRKAnQGqeM51SnQTaQ1pu2DQwJcedrkoCCZLuTqceTaI3mEBmtjtdmUhvQjJTs9zpigQyE5OZljWVvMwU8GFzqelpiVy+ahYzM1O8DsXEMRluhAQR2a6qq4ebNxGKi4t169atE71b4yMisk1Vi8dxex8ALlbVG93pa4HTVfVTYev8BfiWqj7rTj8B3Kyq2wbaJlhZN2M33mU9Wqysm7EarKxH8mC8RUTWhm3oVKBtPIMzJo6VA4Vh0wU47ddGuo4xxphxFMkjqs8AvxeRvgvyTOCDUYvImPiyBVgoInOBI8BVwIf6rfMg8Em3fc46oGGo9jfGGGPGbthHVAAikgAsBgTYq6pd0Q5skDhqgEODLM4C/NbvvB+PCbw9rjmqmj2eGxSRy4Af4rwmfqeq/peIfAxAVW91XxP/KXAJzmviN6jqkHXyQ5R1KxPxxVdlPRqsrPtGzJX1SNrgJAP/DKwHFHgGuDXWBtwUka3x8Lx5JPx4TODf45oIfv3b2XGZ/vz6t7PjmjiRPKL6FU7fNz9xp68G7gZiargGY4wxxpg+kSQ4i1V1Vdj0UyJSEq2AjDHGGGPGKpK3qF4RkTP6JkRkHfBc9EIatdu8DiAK/HhM4N/jmgh+/dvZcZn+/Pq3s+OaIJG0wdmD08D4sDtrNrAH6AVUVVdGNUJjjDHGmBGKJMGZM9RyVR3srSZjjDHGGE9E9Jq4McYYY0w8ifshfkXkEhF5TUT2iciXvI5ntESkUESeEpE9IrJLRD7tzp8uIptE5A335zSvYx0pEQmKyCsi8pA7HffH5AUr67HPyvr4sLIe++KhrMd1guOObH4LcCmwFLhaRJZ6G9WodQOfU9UlwBnAJ9xj+RLwhKouBJ5wp+PNp3HabfXxwzFNKCvrccPK+hhZWY8bMV/W4zrBAU4H9qnqAVXtBO4FrvA4plFR1UpVfdn93oRTcPJxjueX7mq/BN7tSYCjJCIFwDuAn4fNjutj8oiV9RhnZX3cWFmPcfFS1uM9wckHysKmy915cU1EioA1wItAbt+4Re7PHA9DG40fAjfjvHXXJ96PyQtW1mPfD7GyPh6srMe+HxIHZT3eExwZYF5ct5oWkSnAH4DPqGqj1/GMhYi8E6hW1W1ex+IDVtZjmJX1cWVlPYbFU1mPpCfjWFYOFIZNFwAVg6wb89xBTf8A3KOqf3RnV4nITFWtFJGZQLV3EY7Y2cDl7mCUyUCGiPya+D4mr1hZj21W1sePlfXYFjdlPd5rcLYAC0VkrogkAlcBD3oc06iIiAB3AHtU9fthix4ErnO/Xwc8MNGxjZaqfllVC1S1COff5klV/TBxfEwesrIew6ysjysr6zEsnsp6XNfgqGq3iHwSeBQIAneq6i6Pwxqts4FrgR0ist2d9xXg28BGEfkITm/Sfhjk1I/HFFVW1uOWH48pqqysx62YOybr6M8YY4wxvhPvj6iMMcYYY05iCY4xxhhjfMcSHGOMMcb4jiU4xhhjjPEdS3CMMcYY4zuW4BhjjDHGdyzBMcYYY4zvWIJjjDHGGN+xBMcYY4wxvmMJjjHGGGN8xxIcY4wxxviOJTjGGGOM8R1LcGKciKiILPA6DmMmAxH5ioj83Os4jDFjZwnOKIlIqYi0iUhz2OenXsflNRH5hoj82us4zPgRkQ+JyFa3jFeKyF9FZL3XcY2ViGwQkfLwear636p6o1cxGf8RkS+LyMP95r0xyLyrJjY6f7MEZ2zepapTwj6f9DogY8aTiHwW+CHw30AuMBv4X+AKD8MyJp5sBs4WkSCAiOQBCcDafvMWuOuacWIJzjgTketF5DkR+YGI1IvIARE5y51fJiLVInJd2Pq/EJFbRWSTiDSJyNMiMmeQbWeKyK9EpEZEDonI10QkICJJInJMRFaErZvj1jBl992pisjN7v4rReTdInKZiLzu/u5Xwn43ICJfEpH9IlInIhtFZLq7rMh9bHadiBwWkVoR+aq77BLgK8AH3bv9kmj9nU30iUgm8J/AJ1T1j6raoqpdqvpnVf2CW+5+KCIV7ueHIpLk/m5fmftcWJm7IWzbl4nIbrfMHxGRz7vzrxeRZ/vFceIxrXu+/K9bi9Tsnmt57r6Pi8heEVkT9rul7h30bnf5XSKSLCJpwF+BWWE1sLP610CKyOUisss9l/8mIkv6bfvzIvKqiDSIyO9EJDk6/xomjm3BSWhWu9PnAE8Br/Wbtx+4WET2uOfFARH5p/ANudfwSvd8u7HfuZEkIt9zr8tV7v8rKRNwfDHLEpzoWAe8CswAfgPcC5yGk6F/GPipiEwJW/8a4P8BWcB24J5BtvsTIBOYB5wL/ANwg6p2uPv4cNi6VwOPq2qNO50HJAP5wL8Dt7vrnwq8Dfh3EZnnrvsvwLvdfcwCjgO39ItlPbAYuMD93SWq+gjOnf7v3BqtVUP9kUzMOxOnzNw/yPKvAmfgXKRXAacDXwtbnodTXvOBjwC3iMg0d9kdwD+pajqwHHhyBHFd6e4nC+gAngdedqfvA77fb/1rgIuB+cAi4Guq2gJcClSE1cBWhP+SiCwCfgt8BsgGHgb+LCKJ/WK5BJgLrASuH8FxmElAVTuBF3GSGNyfzwDP9pu3GagG3glkADcAPxCRtXDiBvKzwIU4/5ec229X/4NTvle7y/uu9ZOXqtpnFB+gFGgG6sM+H8W5wL0Rtt4KQIHcsHl1wGr3+y+Ae8OWTQF6gEJ3WnEKaxDnYr40bN1/Av7mfl8HlAEBd3orcKX7fQPQBgTd6XR3u+vCtrUNeLf7fQ9wQdiymUAXEAKK3N8tCFv+EnCV+/0bwK+9/vexz7iU8WuAo0Ms3w9cFjZ9MVDqfu8rc6Gw5dXAGe73w275zei3zeuBZ/vNU2CB+/0XwO1hyz4F7AmbXgHUh02XAh8Lm74M2B8WY3m/fZ0ov8C/ARvDlgWAI8CGsG1/OGz5d4Bbvf53s0/sfdxydb/7vQRYiJMYh8+7boDf+xPwaff7ncC3wpYt4M3/HwRoAeaHLT8TOOj1sXv5sRqcsXm3qk4N+9zuzq8KW6cNQFX7zwuvwSnr+6KqzcAxnJqTcFlAInAobN4hnCwdVX0Rp4CfKyKn4BT6B8PWrVPVnvCYBoizL6Y5wP1utXw9TsLTg9MGo8/RsO+t/Y7H+EMdkCUioUGWz+Lk8hhebutUtTtsOrycvA8n2TjkPpY9cwRx9S+3Q51bEHZ+DRDjUN5yfKra624rP2wdOw9MJDYD690azGxVfQP4O3CWO285sFlELhWRF9xmA/U450iWu41ZvLUsh3/PBlKBbWHX7Ufc+ZOWJTixobDvi/voajpQ0W+dWpxalPD2ObNx7ij7/BLnsdO1wH2q2j7KeMqAS/slb8mqemTY33TuKIw/PA+04zyuHEgFJ5fH/uV2QKq6RVWvAHJw7lI3uotacC7UwInGl2NVGPY9PMbhyupbjk9ExN1WJOeBMeGex3lcexPwHICqNuKUsZvcnxXAH4Dv4dT4T8V5LCruNiqBgrBthpfrWpzkflnYNTtTVSd1wm0JTmy4TETWu8/2/x/woqqGZ+e4tS8bgf8SkXRxGiJ/Fgh/Jftu4D04Sc6vxhDPre5+5gCI01A50rdmqoAiEbGyFedUtQHnGf4tbqP0VBFJcO8yv4PTPuVrbvnIctcdtosAEUkUkWtEJFNVu4BGnBpCcKrql4nIarfB7jfG4VA+ISIF4jSU/wrwO3d+FTDDbUw9kI3AO0TkAhFJAD6H85j47+MQk5lEVLUNp9nAZ3Ha3/R51p23GaeGPgmoAbpF5FLg7WHrbgRuEJElIpJKWPsat3bxdpw2OzkAIpIvIhdH76hin/0nNDZ/lrf2gzNYY8zh/Ab4Os6jqVNx2j4M5FM4d7gHcE6M3+A8lwVAVctxGlsqbz2JRupHOI+3HhORJuAFnDY+kfi9+7NORF4eQwwmBqjq93EuwF/DufCWAZ/EqXX5Js5F+1VgB07Z+2aEm74WKBWRRuBjuA3kVfV1nDe3HgfewCnnY/Ub4DGc8+ZAX4yquhcnSTvgVuu/5dGVqr7mxvUTnDvkd+F0DdE5DjGZyedpnBrL8DL9jDtvs6o24bzgsRHnxY4PEdbMQFX/CvwY5w2sfTi1QuAk3QBfdOe/4J5Xj+O8CDJpidsYyXhERH6B09Dxa8OtG+H27sR5M2RctmdMPBORUuBGVX3c61iMGU9ulwU7gaR+bd2My2pwfEREioD34ryCa4wxxkdE5D3uI95pOK+F/9mSm8FZguMTIvL/cLL576rqQa/jMcYYM+7+CedR8X6cdmsf9zac2GaPqIwxxhjjO1aDY4wxxhjfGawDr5iUlZWlRUVFXodh4ti2bdtqVTXmO7+ysm7Gysq6mSwGK+ueJjju2Bo/whmG4Oeq+u2h1i8qKmLr1q0TEpvxJxE5NPxaUdmvlXUzoaysm8lisLLu2SMqcYaJvwVnwLulwNUistSreIyJFivrZrKwsm5iiZc1OKcD+1T1AICI3AtcAewe7Bfq6urYvn07q1evpqenh7vvvpu1a9eycuVKurq6uOeeeyguLmb58uW0t7dz7733sm7dOpYsWUJraysbN27kzDPPZPHixTQ3N3Pfffexfv16FixYQENDA/fffz/nnHMO8+bN4/jx4zzwwANs2LCBoqIiamtreeihh7jgggsoLCykurqahx9+mIsuuoj8/HyOHj3KI488wiWXXEJeXh7l5eU8tmkT68+7iOT06ZSVlfHKi8+w5PQNJE2ZSk1lOaU7XmLOqecSSkmnoaqco6+9Qv7qcwglT6Gpuoy6/TvIWXUOwaRUWqoOU39oF9krziWYmExr9SEaD+8le/V5SDCR1qqDNJe/zozVFxAIhmipPEBLxRtkr7kICQRoqdhHS+V+ck51OrZsPvI6bVWlZK25yJkuf4322nKyVl/gTJftoeP4UWasPA+ApsO76GyoZcYKZwDbptKddDUfY/pyZzDcxoOv0t3ayPRl653pA9vpaW9l2tKzAGjY/zK9XR1MO8UZcqjhjW1obzdTFzv9B9a/vgWAqYtOc6ZfexEJhMhceCoAx/c+TyAhicz5a53p3X8nmJxKxrzVABzb9Syh1Awy5q4kLSlEwuGXyMvLY/16J56NGzdSUBDey/mEmhRl/ezzLiI5fRplZWVsf/HZQct6fVUZVa9tZ9bqc0iwsj6msj4tNRE9+IKV9Qku62dtuIiUjIHL+sEdL1E0QFkPJafRXF1uZX2UZT03I5mON/4+orLuZYKTz1sHCytngN5yReQmnLE6yM/P77/YM13dvbR29vDsvlqO7mmlsrISjjbx4D3bONKRhLQdZ7Uc48evv0C9ppATaGZtqIHbNpbQqMnkBppYE2rijgd20axJzAw0sirUzC8e3kuLJpIfaGBFqJm7/voabSRQGKhnWaiFux59jQ4SmBM4zpJQC3f8dS9dhJgbPMbiYAu3/3UPPQSZF6xjUbCF2x7ejRJgQbCWBcEWfvaQc51ZFKxhbrCFW/+yB4BTgtUUBptPTC8NVjEz2MQT7vTyUDXZgWaeOuxMrwhVM11aefqQM70qVEOmtLO51JleE6olTTp59qAzfWqojiTp5u/7nenTQnUEpZcX9jnTpyccA+ClN5zpMxKO06MBtrzuTJ+VUE+Hhti215len1BPi7byyh5n+pyEBhq0g5Lde5iZmcznFkbln3204rqsd3T10NLZzeY3aqna08LRykqobOSBXztlPdB2nNUBK+telPVFuVO4KXxEIu/Ff1nv6Gbz6zVU7m7h6NFKpLKRB369lfKOZIJtxyIq63daWQfGt6yfOmcaV4+wRZlnr4mLyAeAi1X1Rnf6WuB0Vf3UYL9TXFysXj2r7e1VXik7zuN7qtlaeoyS8gY6u3tPLM9OTyIvI5ncjCSmpyUyLTWRjJQEpiSFSE0MkpoYIjkhQFIoSGIoQEJQSAgGCAWFUEAIBgIEBAIiBAJCQEBwfuLOF0Dcn33EnXjL3PAVhiEjWDdeCJCenDDwMpFtqlo8ofHEWVnv6VW2lB7jqb3VbCk9xs4jjXT2OGVdBLKnJDEzM5mcjGSmpyYyNS2BzBNlPURKQpDkhADJCX1lPUAoIG8p60ERRLCyPkYBEaYkDXyfamV9eN09vbxw4BhPvVbN1kPH2XWkge5e5//EgEBOejJ5mcnkpDvX9ampiWSkhJiSFCItMURKYrDfdd0p6wnBAMEAVtbHUVCEtBGWdS9rcMp562ioBUQ4EvFEqmxo467nSvlzSQWVDe2EAsKy/EyuPWMOKwsyWZAzhXlZU0hJDHodqoldcVHWy461csezB/nLjkpqmjpIDAZYUZDJ9WcXsSI/k4W5UyiakUZygpV1M6i4KOsHapq587mDPLzjKMdaOkkKBVhVMJWPnjOPFfnOdX3OjFSSQlbW45mXCc4WYKGIzAWOAFfhDC4WEyrq27jlqX38fms5vapsWJzNzZcs5sIluYPWDhgziJgu64frWvnJk2/wx1eOEBThgiU5vGPlTM5bnDPoHZMxg4jpsr6vupmfPPkGfy6pICEY4O3L8njHiplsWJxtibsPeXb1UtVuEfkk8CjO64R3quour+IJ96dXjvBvf9pJe3cPHygu5OPnzqdweqrXYZk4FatlXVX57Utl/OdDu1CFa8+Yw8fOnU9eZrLXoZk4FatlvbdXufO5g3znkdcIBYWPnjOPG9fPIzs9yevQTBR5enumqg8DD3sZQ7jmjm6+dv8O/rS9guI50/jBB1dbYmPGRayV9YbWLm7+QwmP7qpi/YIsvveBVZbYmHERa2W9trmDz24sYfPrNVy4JJdvv28FWVMssZkMrP7Z1drZzT/etYVth4/zrxcu4hPnzScUtJEsjP80tnfx4TteZO/RRr562RI+sn4ugYAPWyWaSa+uuYOrbnuBsmOtfPPdy7lm3WzEjy1wzYAswQHau3r4p7u3sfXQMX589RreuXKW1yEZExUtHd3ccNcW9h5t5LZriznvlByvQzImKhpau7j2jpcoP97KL//xdM6YN8PrkMwEm/QJTk+v8snfvMIzb9TyvQ+ssuTG+FZndy8f/dVWtpfVc8uH1lhyY3yrrbOH6+56iX3Vzfz8umJLbiapSZ/g3PXcQR7fU8V/XL6M95/qWe+fxkTdLU/t4+/76/j+lau4ZPlMr8MxJmq+++hrbC+r52fXnso5i2J+vFETJZO6kUlpbQvfe+w1LlySwz+cOcfrcIyJmj2Vjdzy1D7esyaf9661RN7417ZDx7jr7wf5hzPncPGyPK/DMR6atAlOb6/yxT+8SkIwwDffvcIanhnf6u7p5Qv3lTA1NYF/f6eNe2j8q72rhy/c9yqzMlO4+ZJTvA7HeGzSJjj3vHSYFw8e42vvWGKvxxpfu+2ZA+w80sh/XrGcaWmJXodjTNT86Ik3OFDTwrfeu2LQISzM5DEpE5z2rh5+sOl1zpw3gyuLY2ukOmPGU0NbF7c8uY+3L83lshXW7sb4V3VjO3c8c5D3rS2wdjcGmKQJzh9eLudYSyefvnChPZoyvvbblw7T0tnDpy+MreHVjRlvv3y+lK7eXj51/gKvQzExYtIlOL29yh3PHGRFfibr5k73Ohxjoqazu5dfPFfKWfNnsGxWptfhGBM1rZ3d/PqFw7x9aS5FWWleh2NixKRLcJ7cW82B2hZufNtcq70xvvaXHRUcbWzno2+b53UoxkTVfdvKaWjrsrJu3mLSJTi3P3OAWZnJ1h7B+JqqcvvmgyzImcK51h7B+FhPr3LHswdZXTiVU+dM8zocE0MmVYKzo7yBFw8e44az55Jg40wZH3t+fx27Kxu50caZMj63aXcVh+pa+ejb5lmtvHmLSfW//O+3lZGcEOCDp9ubU8bffr+tnKmpCbx7Tb7XoRgTVfdtK2NmZjIXL8v1OhQTYyZNgqOqPL67inMWZpORnOB1OMZETVdPL0/ureaCU3JJTgh6HY4xUdPa2c0zb9Ry8bI8QlYrb/qZNCViV0UjFQ3tXLTUsnzjb1tKj9HQ1mVl3fjeM2/U0tHdy9utrJsBTJoE57HdVQQEzrcRlI3PbdpdRWIowNsWZnkdijFRtWl3FRnJIU6zLj/MACZNgrNpdxWnzpnGjClJXodiJhER+a6I7BWRV0XkfhGZGs39qSqbdlexfkEWadZVvfGxnl7lyb3VnHdKjr00YgY0KUpF2bFW9lQ2WpW98cImYLmqrgReB74czZ3tPdpE+fE2K+vG97YdOs6xlk4r62ZQkyLBeWJPFQAXLc3zOBIz2ajqY6ra7U6+ABREc3+P73bK+gVL7FGs8bfH91SREBTr58kMalIkOJv2VDE/O4251oW38dY/An8dbKGI3CQiW0Vka01Nzah2sGlPFasLp5KTnjzaGI2JeX2PYs+YN4N0eyvWDML3CU5DWxcvHjhmtTcmakTkcRHZOcDnirB1vgp0A/cMth1VvU1Vi1W1ODt75HelRxvaebW8warsje/tr2nhYG2LvT1lhuT7VoivHD5Od69yziJ7o8REh6peONRyEbkOeCdwgapqtOLYUnoMwKrsje+9dNAp6+dYWTdD8KQGZyLfLCkpa0AEVhZEbRfGDEpELgG+CFyuqq3R3FdJWT1JoQCL89KjuRtjBiQiHxCRXSLSKyLF0dxXSVk901ITmD09NZq7MXHOq0dUE/ZmSUl5PQtzpjDFXpk13vgpkA5sEpHtInJrtHZUUl7P8vxMe2XWeGUn8F5gc7R3VFJez6rCqTb2lBmSJ//rq+pjYZMvAO+P0n4oKau3zv2MZ1R1wUTsp7unlx1HGvjQ6XMmYnfGnERV9wBRTzpaOrp5vaqJi5dZu0oztFi41YvamyXlx9uoa+lkVeHUMYZoTGx7vaqZ9q5eVhVmeh2KMcMay3V955EGehVW23XdDCNqNTgi8jgwUIr9VVV9wF0nojdLgNsAiouLR9RAs6S8HrATwfiflXUzESK5rkdiPK7rKwssmTdDi1qCEwtvlpSU1ZNojS7NJFBSVs9Ua3Rpomy46/pEKClroHB6ig27Y4blSRucsDdLzo3mmyUlZQ0sn5VhjS6N720vq2dVgTW6NP63vayeNbOneh2GiQNe/c8f9TdL+hpdWvsb43etnU6jSyvrxksi8h4RKQfOBP4iIo+O9z5qmjo4Ut9mj2JNRLx6iyrqb5a8Ud1MW1ePnQjG93YeaXQbXVqbBOMdVb0fuD+a+3jVbX9jybyJhG+f3ZSU1QOwyjr4Mz7XV9atM0vjdyVl9QQDwrJZGV6HYuKAfxOc8noyUxKYM8MaXRp/215eT8G0FLKs0aXxue3lDSzKTSc10TpuNcPzb4JT1sDKgkxrdGl879XyequpNL6nqm5Zt0exJjK+THB6e5X9Nc0szrXXw42/tXf1UHasjUVW1o3P1bV0Ut/aZWXdRMyXCc7RxnY6unspykrzOhRjoupQndPLQlGWPYo1/lZa2wLAXLuumwj5MsGxE8FMFgetrJtJoq+s242riZQvE5yDdXYimMmh1Mq6mSRK61oIBoSCaSleh2LihC8TnNLaFpJCAWZmJHsdijFRVVrbwoy0RDKSE7wOxZioKq1tpXBaivVMbyLmy5JSWtfKnBmpBAL2BpXxt9K6Fqu9MZOClXUzUv5McGpbKJphJ4Lxv9LaVivrxvdU1a7rZsR8l+D09iqHjrVao0vje22dPRxtbGeuvUFlfK6muYOWzh67rpsR8V2CU9HQRqe9Im5ijIh8XkRURLLGa5vWwNhMFqW1fd0hWFk3kfNdgnPiRLCqTBMjRKQQuAg4PJ7b7esOwcq68bsTXX9YWTcj4LsEp+8VcavKNDHkB8DNgI7nRq07BDNZHKxrISEozJpqb8aayPkuwSmtbSElIUhuhg08aLwnIpcDR1S1JIJ1bxKRrSKytaamZthtl9a2kJ2exJQkG3jQ+FtpbQuF01MJ2SviZgR8d2UsrW1hzoxUG2TTTBgReRzIG2DRV4GvAG+PZDuqehtwG0BxcfGwtT2lta1WZW8mhYO1LVbWzYj5LsE5WNdig2yaCaWqFw40X0RWAHOBEjfhLgBeFpHTVfXoWPd7sK6F8xZnj3UzxsQ0VeVQXStnLxi39vlmkvBVfV93Ty9lx1qtTYKJCaq6Q1VzVLVIVYuAcmDteCQ3zR3d1DR1WFk3vlfV2EFbV4+VdTNivkpwKurb6epRq8o0vmdvlZjJ4qCVdTNKvnpEZW+VmFjm1uKMC+sDx0wWb5Z169DSjIyvanDe7BfETgTjb31lfY6VdeNzpbUtJAYDzMy0UcTNyPgqwSk71kpKQpDsdHtF3Phb2bE2stOTSE30VSWsMScpO95KwfQUgjZ4shkhXyU4VU0d5GYk2SvixveqmtrJy7BOz4z/VTV2WFk3o+JpgjPe4/NUNbaTk24ngvG/qsYOcqym0sQQEfmuiOwVkVdF5H4RmToe23Wu61bWzch5luBEY3yemqYOcqwHYzMJ1DS1k2N3tSa2bAKWq+pK4HXgy2PdoKpS3dRBrpV1Mwpe1uCM+/g81VaDYyaBrp5e6lo67a7WxBRVfUxVu93JF3A6thyTxrZuOrt7rV2lGRVPEpxojM/T3NFNS2eP1eAY36tt7kAVK+smlv0j8NfBFkZ6Xa9qagew2kozKlF7BWOix+epbnROBBtk0/hddWMHALlWW2km2FDXdVV9wF3nq0A3cM9g24n8ut5X1u26bkYuagnORI/PU93knAj2iMr43Ymybsm8mWCDXdf7iMh1wDuBC1R1zM0Pqq0Gx4zBhHeioao7gJy+aREpBYpVtXYs261ya3CsXYLxuzfLul30TewQkUuALwLnqmrreGyzqrHvxtWu62bkfNMPTs2Ju1q76Bt/q27qQASypiR6HYox4X4KpAObRGS7iNw61g1WN7UzJSlEWpJ1aGlGzvNSM17j81Q3dZAUCpCR7PkhGRNVNU3tzEhLIhT0zf2J8QFVXTDe26xusv6ezOj55gpZ1dhOjvVibCYB6+TPTBbVje32irgZNd8kONWNHfZWiZkUqpva7W1BMylYJ39mLPyT4DS121slZlKobuywBsbG91TVLet2XTej458Exy76ZhLo6VVqm21IEuN/TR3dtHVZ561m9HyR4LR19tDU0W0ngvG9uuYOetXeFjT+d6KTPyvrZpR8keCc6AzKanBMDBKRT4nIayKyS0S+M5ZtvdmhpSXzxt/6ruvWyNiMli/eqbbOoEysEpHzgCuAlaraISI5w/3OUKxDSzNZVDda7/RmbHxVg2NVmSYGfRz4tqp2AKhq9Vg21leDY2Xd+N2b13VL5s3o+CPBsRocE7sWAW8TkRdF5GkROW2wFSMZYbmvrGdNsbJu/K26sYOUhCBTrBdjM0q+KDlVTe0kBgNMTU3wOhQzCQ01wjLOOTYNOAM4DdgoIvMGGogwkhGWq5ramZ6WSGLIF/cmxgyqqqnDOm81Y+KLBKemsYPsdDsRjDeGGmFZRD4O/NFNaF4SkV4gCxi4imYY1i+ImSyqG9ut81YzJr64Daxusn5BTMz6E3A+gIgsAhKB2tFurKap3V4RN5NCTVMH2XZdN2PgiwSnqrHd7mpNrLoTmCciO4F7gesGejwVKRuHykwWdl03Y+WLR1TVTR2cOX+G12EYcxJV7QQ+PB7b6nV7Mba3SozftXR009LZY28LmjGJ+xqc9q4eGtq6LNM3vnestZPuXrV+QYzvWYeWZjzEfYJT02SdQZnJwTr5M5PFm2Xdrutm9OI+wTkxTINV2xufO3FXa9X2xufe7NDSrutm9OK+Dc7ivAx+d9MZnJKX4XUoxkTV2tnTuPemM1icl+51KMZE1dsWZPGbj66jcHqq16GYOBb3Cc6UpBDr5lkDY+N/mSkJnGFl3UwC09ISOWt+ltdhmDgX94+ojDHGGGP6swTHGGOMMb4jY+hzbMKJSA1waJDFWYyhh9gY5cdjAm+Pa46qZnu074gNUdatTMQXK+vDsLLuGzFX1uMqwRmKiGxV1WKv4xhPfjwm8O9xTQS//u3suEx/fv3b2XFNHHtEZYwxxhjfsQTHGGOMMb7jpwTnNq8DiAI/HhP497gmgl//dnZcpj+//u3suCaIb9rgGGOMMcb08VMNjjHGGGMMYAmOMcYYY3wo7hMcEblERF4TkX0i8iWv4xktESkUkadEZI+I7BKRT7vzp4vIJhF5w/05zetYR0pEgiLyiog85E7H/TF5wcp67LOyPj6srMe+eCjrcZ3giEgQuAW4FFgKXC0iS72NatS6gc+p6hLgDOAT7rF8CXhCVRcCT7jT8ebTwJ6waT8c04Sysh43rKyPkZX1uBHzZT2uExzgdGCfqh5Q1U7gXuAKj2MaFVWtVNWX3e9NOAUnH+d4fumu9kvg3Z4EOEoiUgC8A/h52Oy4PiaPWFmPcVbWx42V9RgXL2U93hOcfKAsbLrcnRfXRKQIWAO8COSqaiU4JwuQ42Foo/FD4GagN2xevB+TF6ysx74fYmV9PFhZj30/JA7KerwnODLAvLh+711EpgB/AD6jqo1exzMWIvJOoFpVt3kdiw9YWY9hVtbHlZX1GBZPZT3kdQBjVA4Uhk0XABUexTJmIpKAcxLco6p/dGdXichMVa0UkZlAtXcRjtjZwOUichmQDGSIyK+J72PyipX12GZlffxYWY9tcVPW470GZwuwUETmikgicBXwoMcxjYqICHAHsEdVvx+26EHgOvf7dcADEx3baKnql1W1QFWLcP5tnlTVDxPHx+QhK+sxzMr6uLKyHsPiqazHdQ2OqnaLyCeBR4EgcKeq7vI4rNE6G7gW2CEi2915XwG+DWwUkY8Ah4EPeBPeuPLjMUWVlfW45cdjiior63Er5o7JhmowxhhjjO/E+yMqY4wxxpiTWIJjjDHGGN+xBMcYY4wxvmMJjjHGGGN8xxIcY4wxxviOJTjGGGOM8R1LcIwxxhjjO/8/SOdrPMKcTCQAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 576x288 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plotT = 50\n",
    "tfp = (irf['C'] + ss['C'])/(irf['L'] + ss['L'])/((ss['C']/ss['L']))\n",
    "fig, axes = plt.subplots(2, 3, figsize=(8, 4))\n",
    "ax = axes.flatten()\n",
    "\n",
    "ax[0].plot(100*((irf['N'][0:plotT] + ss['N'])/ss['N']-1))\n",
    "ax[0].axhline(0, color='gray', linestyle=':')\n",
    "ax[0].set_title('Mass of establishments')\n",
    "ax[0].set_ylabel('pct deviation from ss')\n",
    "\n",
    "ax[1].plot(100*((irf['mu'][0:plotT] + ss['mu'])/ss['mu']-1))\n",
    "ax[1].axhline(0, color='gray', linestyle=':')\n",
    "ax[1].set_title('Markup')\n",
    "\n",
    "ax[2].plot(100*((tfp[0:plotT] - 1)))\n",
    "ax[2].axhline(0, color='gray', linestyle=':')\n",
    "ax[2].set_title('TFP')\n",
    "\n",
    "ax[3].plot(100*((irf['L'][0:plotT] + ss['L'])/ss['L']-1))\n",
    "ax[3].axhline(0, color='gray', linestyle=':')\n",
    "ax[3].set_title('Employment')\n",
    "\n",
    "ax[4].plot(100*((irf['C'][0:plotT] + ss['C'])/ss['C']-1))\n",
    "ax[4].axhline(0, color='gray', linestyle=':')\n",
    "ax[4].set_title('Consumption')\n",
    "\n",
    "ax[5].plot(100*((irf['w'][0:plotT] + ss['w'])/ss['w']-1))\n",
    "ax[5].axhline(0, color='gray', linestyle=':')\n",
    "ax[5].set_title('Wage')\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "c87a2dd7",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = [100*((irf['N'][0:plotT] + ss['N'])/ss['N']-1),\n",
    "       100*((irf['mu'][0:plotT] + ss['mu'])/ss['mu']-1),\n",
    "       100*((tfp[0:plotT] - 1)),\n",
    "       100*((irf['L'][0:plotT] + ss['L'])/ss['L']-1),\n",
    "       100*((irf['C'][0:plotT] + ss['C'])/ss['C']-1),\n",
    "       100*((irf['w'][0:plotT] + ss['w'])/ss['w']-1)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "df1115a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = np.array(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "a2fe712a",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt(\"translog_irf.csv\", np.transpose(data), delimiter=\",\", fmt = \"%g\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
